


% Load Matrices and Estimates    
load('Mats_T2C.mat')
if Net_OUT == 0
    load(['PE_Estimates_',num2str(sim_tol),'_OR.mat']);
else
    load(['PE_Estimates_',num2str(sim_tol),'_OUT.mat']);
end    

Params_Together_b = [];
V_Together_b = [];

SIM_SUM = size(PE_Chain_Out.PE_index,2)

for z=1:SIM_SUM
    z
    
    y_z = struct('y1_educ_m1',PE_Chain_Out.y1_educ_m1_b(:,z),'y1_educ_m2',PE_Chain_Out.y1_educ_m2_b(:,z),'y1_educ_m3',PE_Chain_Out.y1_educ_m3_b(:,z),'y1_educ_m4',PE_Chain_Out.y1_educ_m4_b(:,z),...
        'y1_gender_m1',PE_Chain_Out.y1_gender_m1_b(:,z),'y1_gender_m2',PE_Chain_Out.y1_gender_m2_b(:,z),'y1_gender_m3',PE_Chain_Out.y1_gender_m3_b(:,z),'y1_gender_m4',PE_Chain_Out.y1_gender_m4_b(:,z));

    Params_z = [PE_Chain_Out.BGD_hat_b(:,z); PE_Chain_Out.A_hat_b(:,z); ...
        PE_Chain_Out.alpha_educ_m1_b(:,z); PE_Chain_Out.alpha_educ_m2_b(:,z); PE_Chain_Out.alpha_educ_m3_b(:,z); PE_Chain_Out.alpha_educ_m4_b(:,z); ...
        PE_Chain_Out.alpha_gender_m1_b(:,z); PE_Chain_Out.alpha_gender_m2_b(:,z); PE_Chain_Out.alpha_gender_m3_b(:,z); PE_Chain_Out.alpha_gender_m4_b(:,z)];        

    [A, B, Var_Together_z] = GMM_Fn_Mar2022(Params_z,Cons_T2C,PE_Chain_Out.G_ij_b(:,z),y_z);
    Params_Together_b = cat(3,Params_Together_b,Params_z);
    V_Together_b = cat(3,V_Together_b,Var_Together_z);

end



dim_BGD = 4*Cons_T2C.dimX+3;

Qbar_Together = zeros(dim_BGD+Cons_T2C.n+2*(60 ),1);
Ubar_Together = zeros(dim_BGD+Cons_T2C.n+2*(60 ),dim_BGD+Cons_T2C.n+2*(60 ));

for z=1:SIM_SUM
    Qbar_Together = Qbar_Together + Params_Together_b(:,:,z)/SIM_SUM;
    Ubar_Together = Ubar_Together + V_Together_b(:,:,z)/SIM_SUM;
end

Bbar_Together = zeros(dim_BGD+Cons_T2C.n+2*(60 ),dim_BGD+Cons_T2C.n+2*(60 ));
for imp=z:SIM_SUM
    Bbar_Together = Bbar_Together + (Params_Together_b(:,:,z) - Qbar_Together)*(Params_Together_b(:,:,z) - Qbar_Together)'/(SIM_SUM-1);
end


VarEst.Var_Together = Ubar_Together + (1 + 1/SIM_SUM)*Bbar_Together;



%%% Gather all final stuff together
% Point estimates
VarEst.BGD_hat = Qbar_Together(1:dim_BGD,:);
VarEst.A_hat = Qbar_Together(dim_BGD+1:dim_BGD+Cons_T2C.n,:);

VarEst.alpha_educ_m1 = Qbar_Together(dim_BGD+Cons_T2C.n+1:dim_BGD+Cons_T2C.n+9,:);
VarEst.alpha_educ_m2 = Qbar_Together(dim_BGD+Cons_T2C.n+10:dim_BGD+Cons_T2C.n+24,:);
VarEst.alpha_educ_m3 = Qbar_Together(dim_BGD+Cons_T2C.n+25:dim_BGD+Cons_T2C.n+39,:);
VarEst.alpha_educ_m4 = Qbar_Together(dim_BGD+Cons_T2C.n+40:dim_BGD+Cons_T2C.n+60,:);



VarEst.alpha_gender_m1 = Qbar_Together(dim_BGD+Cons_T2C.n+60 +1:dim_BGD+60 +Cons_T2C.n+9,:);
VarEst.alpha_gender_m2 = Qbar_Together(dim_BGD+Cons_T2C.n+60 +10:dim_BGD+60 +Cons_T2C.n+24,:);
VarEst.alpha_gender_m3 = Qbar_Together(dim_BGD+Cons_T2C.n+60 +25:dim_BGD+60 +Cons_T2C.n+39,:);
VarEst.alpha_gender_m4 = Qbar_Together(dim_BGD+Cons_T2C.n+60 +40:dim_BGD+60 +Cons_T2C.n+60,:);




% Variances
VarEst.Var_BGD_hat = VarEst.Var_Together(1:dim_BGD,1:dim_BGD);
VarEst.Var_A_hat = VarEst.Var_Together(dim_BGD+1:dim_BGD+Cons_T2C.n,dim_BGD+1:dim_BGD+Cons_T2C.n);

VarEst.Var_alpha_educ_m1 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+1:dim_BGD+Cons_T2C.n+9,dim_BGD+Cons_T2C.n+1:dim_BGD+Cons_T2C.n+9);
VarEst.Var_alpha_educ_m2 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+10:dim_BGD+Cons_T2C.n+24,dim_BGD+Cons_T2C.n+10:dim_BGD+Cons_T2C.n+24);
VarEst.Var_alpha_educ_m3 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+25:dim_BGD+Cons_T2C.n+39,dim_BGD+Cons_T2C.n+25:dim_BGD+Cons_T2C.n+39);
VarEst.Var_alpha_educ_m4 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+40:dim_BGD+Cons_T2C.n+60,dim_BGD+Cons_T2C.n+40:dim_BGD+Cons_T2C.n+60);

VarEst.Var_alpha_gender_m1 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+60 +1:dim_BGD+60 +Cons_T2C.n+9,dim_BGD+Cons_T2C.n+60 +1:dim_BGD+60 +Cons_T2C.n+9);
VarEst.Var_alpha_gender_m2 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+60 +10:dim_BGD+60 +Cons_T2C.n+24,dim_BGD+Cons_T2C.n+60 +10:dim_BGD+60 +Cons_T2C.n+24);
VarEst.Var_alpha_gender_m3 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+60 +25:dim_BGD+60 +Cons_T2C.n+39,dim_BGD+Cons_T2C.n+60 +25:dim_BGD+60 +Cons_T2C.n+39);
VarEst.Var_alpha_gender_m4 = VarEst.Var_Together(dim_BGD+Cons_T2C.n+60 +40:dim_BGD+60 +Cons_T2C.n+60,dim_BGD+Cons_T2C.n+60 +40:dim_BGD+60 +Cons_T2C.n+60);




VarEst.chains = chains;
VarEst.SIM_SUM = SIM_SUM;
VarEst.PE_sim_reps = PE_Chain_Out.PE_sim_reps;
VarEst.b_gap = PE_Chain_Out.b_gap
VarEst.burnin = PE_Chain_Out.burnin

if Net_OUT == 0
    save(['Var_Estimates_',num2str(sim_tol),'_OR.mat'],'VarEst')
else
    save(['Var_Estimates_',num2str(sim_tol),'_OUT.mat'],'VarEst')
end
    

    
